#####REPLICATION FILE FOR "ETHNIC BARRIERS TO CIVIL RESISTANCE"#####

##Packages
library(dplyr)
library(multiwayvcov)   #v. 1.2.2
library(lmtest)
library(car)
library(nnet)
library(simcf)
library(MASS)
library(tile)
library(RColorBrewer)
library(stargazer)
library(zeligverse)

#The best way to install "simcf" and "tile" is to download directly from Chris Adolph's Github:
#install.packages("devtools")
#library(devtools)
#install_github("chrisadolph/simcf")
#install_github("chrisadolph/tile")

###Load Data###
EGC<-read.csv("EthnicGroupsinContention.csv", na.strings=c("NA", ""))
EBCR<-read.csv("EBCR_Replication_Data.csv", na.strings=c("NA", ""))

###Campaign-Level Descriptive Statistics###(Table 1)
MERGE<-merge(EGC, EBCR, by=c("YEAR","EPR_ID"), all.x=TRUE)
MERGE<-MERGE[!is.na(MERGE$EPR_ID),] #Drop Missing, 382 cases remain
CAMP_LEVEL<-MERGE%>%
  group_by(NAVCO2ID, CAMPAIGN, YEAR, LCCODE, LOCATION, TCCODE, TARGET, CAMP_GOALS, PRIM_METHOD, NAVCO1, CDIVERS_ETHNIC, MOVEMENT)%>%
  summarise(TOT_GROUPS=sum(PARTICIPATION.x), 
            TOT_GROUP_SIZE=sum(EPR_GROUPSIZE),
            ANY_EGIP=max(EPR_STATUS_EGIP),
            ANY_SRPL=max(EPR_STATUS_SRPL),
            ANY_JP=max(EPR_STATUS_JP),
            ANY_EXCL=max(EPR_STATUS_EXCL),
            ANY_IRR=max(EPR_STATUS_IRR),
            ANY_CLAIM=max(CLAIM.x),
            INIT_GROUP_SIZE=sum(EPR_GROUPSIZE[INITIATION.x==1]),
            INIT_EGIP=sum(EPR_STATUS_EGIP[INITIATION.x==1]),
            INIT_SRPL=sum(EPR_STATUS_SRPL[INITIATION.x==1]),
            INIT_JP=sum(EPR_STATUS_JP[INITIATION.x==1]),
            INIT_EXCL=sum(EPR_STATUS_EXCL[INITIATION.x==1]),
            INIT_IRR=sum(EPR_STATUS_IRR[INITIATION.x==1]),
            INIT_CLAIM=sum(CLAIM.x[INITIATION.x==1]))%>%
  mutate(MULTI_GROUP=ifelse(TOT_GROUPS>1,1,0))

CAMP_STATS<-CAMP_LEVEL%>%
  mutate(STRATEGY=ifelse(PRIM_METHOD==1, "NONVIOLENT", "VIOLENT"))%>%
  group_by(STRATEGY)%>%
  summarise(
    OBSERVATIONS=n(),
    NAVCO_DIVERSE=mean(CDIVERS_ETHNIC, na.rm=TRUE),
    MULTI_GROUP=mean(MULTI_GROUP, na.rm=TRUE),
    AVG_GROUPS=mean(TOT_GROUPS),
    AVG_SIZE=mean(TOT_GROUP_SIZE, na.rm=TRUE),
    ANY_SRPL=mean(ANY_SRPL),
    ANY_JP=mean(ANY_JP),
    ANY_EXCL=mean(ANY_EXCL),
    ANY_IRR=mean(ANY_IRR),
    ANY_CLAIM=mean(ANY_CLAIM),
    AVG_INIT_SIZE=mean(INIT_GROUP_SIZE, na.rm=TRUE),
    INIT_SRPL=mean(INIT_SRPL),
    INIT_JP=mean(INIT_JP),
    INIT_EXCL=mean(INIT_EXCL),
    INIT_IRR=mean(INIT_IRR),
    INIT_CLAIM=mean(INIT_CLAIM)
  )
CAMP_STATS<-data.frame(t(CAMP_STATS))

###Compare to ACD2EPR### (Table 1)
ACD2EPR<-read.csv(url("https://icr.ethz.ch/data/epr/acd2epr/ACD2EPR-2014.csv"), na.strings=c("NA", ""))
ACD2EPR<-ACD2EPR%>%
  mutate(EPR_ID=gwgroupid, YEAR=from)
ACD_MERGE<-merge(ACD2EPR, EBCR, by=c("YEAR","EPR_ID"), all.x=TRUE)
ACD_MERGE<-ACD_MERGE[!is.na(ACD_MERGE$EPR_GROUP),] #drop missing, 365 cases remain

ACD_CAMP_LEVEL<-ACD_MERGE%>%
  mutate(ACD2EPR_PART=1, ACD2EPR_CLAIM=ifelse(claim>0,1,0))%>%
  group_by(dyadid, YEAR, GW_COUNTRY_ID, statename, sideb)%>%
  summarise(TOT_GROUPS=sum(ACD2EPR_PART), 
            TOT_GROUP_SIZE=sum(EPR_GROUPSIZE),
            ANY_EGIP=max(EPR_STATUS_EGIP),
            ANY_SRPL=max(EPR_STATUS_SRPL),
            ANY_JP=max(EPR_STATUS_JP),
            ANY_EXCL=max(EPR_STATUS_EXCL),
            ANY_IRR=max(EPR_STATUS_IRR),
            ANY_CLAIM=max(ACD2EPR_CLAIM))

ACD_STATS<-ACD_CAMP_LEVEL%>%
  mutate(MULTI_GROUP=ifelse(TOT_GROUPS>1,1,0))%>%
  mutate(STRATEGY= "ACD VIOLENT")%>%
  group_by(STRATEGY)%>%
  summarise(
    OBSERVATIONS=n(),
    MULTI_GROUP=mean(MULTI_GROUP, na.rm=TRUE),
    AVG_GROUPS=mean(TOT_GROUPS),
    AVG_SIZE=mean(TOT_GROUP_SIZE, na.rm=TRUE),
    ANY_SRPL=mean(ANY_SRPL),
    ANY_JP=mean(ANY_JP),
    ANY_EXCL=mean(ANY_EXCL),
    ANY_IRR=mean(ANY_IRR),
    ANY_CLAIM=mean(ANY_CLAIM, na.rm=TRUE))
ACD_STATS<-data.frame(t(ACD_STATS))

#Test for statistically significant differences (Appendix Table 5)
NV_CAMPS<-CAMP_LEVEL[CAMP_LEVEL$PRIM_METHOD==1,]
V_CAMPS<-CAMP_LEVEL[CAMP_LEVEL$PRIM_METHOD==0,]

MULTI_GROUP.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$MULTI_GROUP))
AVG_GROUPS.ttest<-t.test(NV_CAMPS$TOT_GROUPS, V_CAMPS$TOT_GROUPS)
AVG_GROUP_SIZE.ttest<-t.test(NV_CAMPS$TOT_GROUP_SIZE, V_CAMPS$TOT_GROUP_SIZE)
ANY_SRPL.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$ANY_SRPL))
ANY_JP.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$ANY_JP))
ANY_EXCL.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$ANY_EXCL))
ANY_IRR.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$ANY_IRR))
ANY_CLAIM.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$ANY_CLAIM))
AVG_INIT_SIZE.ttest<-t.test(NV_CAMPS$INIT_GROUP_SIZE, V_CAMPS$INIT_GROUP_SIZE)
INIT_SRPL.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$INIT_SRPL))
INIT_JP.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$INIT_JP))
INIT_EXCL.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$INIT_EXCL))
INIT_IRR.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$INIT_IRR))
INIT_CLAIM.chisqtest<-chisq.test(table(CAMP_LEVEL$PRIM_METHOD, CAMP_LEVEL$INIT_CLAIM))


###Group-Level Descriptive Statistics (Table 2) ###


table(EBCR$EPR_STATUS_SRPL, EBCR$INIT_MN_ONSET)
table(EBCR$EPR_STATUS_JP, EBCR$INIT_MN_ONSET)
table(EBCR$EPR_STATUS_EXCL, EBCR$INIT_MN_ONSET)


###Group-Level Regression Analysis###

#Knock out Irrelevant Groups
EBCR<-EBCR[EBCR$EPR_STATUS_IRR==0,]

#Eliminate Venezuelan cases of both violent and nonviolent onset in year
EBCR$INIT_MN_ONSET<-ifelse(EBCR$INIT_MN_ONSET==3, NA, EBCR$INIT_MN_ONSET)
EBCR$MN_ONSET<-ifelse(EBCR$MN_ONSET==3, NA, EBCR$MN_ONSET)
EBCR$RC_MN_ONSET<-ifelse(EBCR$RC_MN_ONSET==3, NA, EBCR$RC_MN_ONSET)

#Stargazer to create summary statistics ofcovariates in Appendix Table 6
stargazer(EBCR)

##Producing Table 3##

#Base Multinomial Model#
BASE.MODEL<-INIT_MN_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG  + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3) 
BASE.DATA<-extractdata(BASE.MODEL, EBCR, na.rm=T)
BASE.MNM<-multinom(BASE.MODEL, data=BASE.DATA, Hess = T, maxit=200)
summary(BASE.MNM)

#Group-level Covariates
GROUP.MNM<-multinom(INIT_MN_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + EPR_TEK_EGIP + EPR_DOWNGRADED5 + HORIZ_INEQ + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3), data=EBCR, maxit=200)
summary(GROUP.MNM)

#Bivariate IV: Exclusion
EXCL.MNM<-multinom(INIT_MN_ONSET ~ EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3), data=EBCR, maxit=200)
summary(EXCL.MNM)


#Stargazer to generate latex code to produce Table 3

stargazer(BASE.MNM, GROUP.MNM, EXCL.MNM,
          title = "Multinomial Logit: Initiating Violent and Nonviolent Campaigns", 
          font.size = "tiny",
          align = TRUE, 
          notes = c("Cubic polynomial for time not displayed.",
                    "GDP and population lagged 1 year."
          ),
          dep.var.caption = "Odd Models: DV = Nonviolent; Even Models: DV = Violent",
          dep.var.labels.include = TRUE,
          covariate.labels = c(
            "EPR Status: Jr. Partner",
            "EPR Status: Excluded",
            "Group Size",
            "Country Population (log)",
            "Country GDP Per Capita (log)",
            "Prior Group NV Participation",
            "Prior Group V Participation",
            "VDEM Polyarchy Index",
            "VDEM Physical Integrity Index",
            "Neighboring Kin in Power",
            "Downgraded",
            "Horizontal Inequality",
            "NV YEARS",
            "NV YEARS2",
            "NV YEARS3",
            "V YEARS",
            "V YEARS2",
            "V YEARS3"
          ),
          colnames = TRUE)


##Table 4 Sensitivity Analysis: Alternative Operationalizations of the DV##

#Regime-Change Campaigns Only
EBCR$INIT_RC_MN_ONSET<-ifelse(EBCR$INIT_MN_ONSET==1 & EBCR$RC_MN_ONSET ==1,1,
                             ifelse(EBCR$INIT_MN_ONSET==2 & EBCR$RC_MN_ONSET ==2,2,0))
RC.MNM<-multinom(INIT_RC_MN_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3), data=EBCR, maxit=200)
summary(RC.MNM)

#Any Participating Group
PART.MNM<-multinom(MN_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL+ EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3), data=EBCR, maxit=200)
summary(PART.MNM)

#Stargazer to generate latex code to produce Table 4

stargazer(RC.MNM, PART.MNM,
          title = "Alternate Dependent Variables", 
          font.size = "tiny",
          align = TRUE, 
          notes = c("Cubic polynomial for time not displayed.",
                    "GDP and population lagged 1 year."
          ),
          dep.var.caption = "Odd Models: DV = Nonviolent; Even Models: DV = Violent",
          dep.var.labels.include = TRUE,
          covariate.labels = c(
            "EPR Status: Jr. Partner",
            "EPR Status: Excluded",
            "Group Size",
            "Country Population (log)",
            "Country GDP Per Capita (log)",
            "Prior Group NV Participation",
            "Prior Group V Participation",
            "VDEM Polyarchy Index",
            "VDEM Physical Integrity Index",
            "NV YEARS",
            "NV YEARS2",
            "NV YEARS3",
            "V YEARS",
            "V YEARS2",
            "V YEARS3"
          ),
          colnames = TRUE)


###Producing Figures 1 and 2###
###Simulating Relative Risks and Expected Values for Figure 1 Using Chris Adolph's simCF package
# Chris Adolph's homemade extractor function for multinom() coef's
coef.multinom <- function(x) {
  nlevel <- length(BASE.MNM$lev)
  ncoef <- length(BASE.MNM$coefnames)
  coef <- x$wts[(ncoef+2):length(x$wts)]
  coef[-((0:(nlevel-2))*(ncoef+1) + 1)]
}
BASE.pe<-coef(BASE.MNM)
BASE.vc <- solve(BASE.MNM$Hess) 

##Simulate parameters from predictive distributions
sims <- 10000
BASE.simbetas <- mvrnorm(sims,BASE.pe,BASE.vc) # draw parameters, using MASS::mvrnorm (matrix with 10000 rows and 2x parameters columns)
BASE.simB <- array(NA, dim = c(sims,16,2))  # re-arrange simulates to array format (like 3 dimensionsal: sims x parameters x outcomes-1)
BASE.simB[,,1] <- BASE.simbetas[,1:16]           #   for MNL simulation (filling in first dimension with parameters for outcome 1)
BASE.simB[,,2] <- BASE.simbetas[,17:32]           # filling in second dimension with parameters for outcome 2)

## Create several comparisons for relative risks moving from 0 to 1 or 1SD below the mean to 1SD above
BASE.FIG1HYP <- cfMake(BASE.MODEL, BASE.DATA, nscen=9)

#Scenario 1: Junior Partner
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "EPR Status: Jr. Partner", scen=1)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "EPR_STATUS_JP",
                         x=1,
                         xpre=0,
                         scen=1)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "EPR_STATUS_EXCL", 
                         x = 0, 
                         xpre=0,
                         scen=1)

# Scenario 2: Excluded
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "EPR Status: Excluded", scen=2)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "EPR_STATUS_EXCL",
                         x=1,
                         xpre=0,
                         scen=2)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "EPR_STATUS_JP",
                         x=0,
                         xpre=0,
                         scen=2)

# Scenario 3: Group size
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "Group Size", scen=3)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "EPR_GROUPSIZE",
                         x=mean(BASE.DATA$EPR_GROUPSIZE) + sd(BASE.DATA$EPR_GROUPSIZE),
                         xpre=mean(BASE.DATA$EPR_GROUPSIZE),
                         scen=3)

# Scenario 4: Population
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "Country Population (log)", scen=4)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "POP_LOG_LAG_EXT",
                         x=mean(BASE.DATA$POP_LOG_LAG_EXT) + sd(BASE.DATA$POP_LOG_LAG_EXT),
                         xpre=mean(BASE.DATA$POP_LOG_LAG_EXT),
                         scen=4)

# Scenario 5: GDP
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "Country GDP per capita (log)", scen=5)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "GDPPC_LOG_LAG_EXT",
                         x=mean(BASE.DATA$GDPPC_LOG_LAG_EXT) + sd(BASE.DATA$GDPPC_LOG_LAG_EXT),
                         xpre=mean(BASE.DATA$GDPPC_LOG_LAG_EXT),
                         scen=5)

# Scenario 6: Prior Group NV Participation
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "Prior Group NV Participation", scen=6)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "PASTNV",
                         x=mean(BASE.DATA$PASTNV) + sd(BASE.DATA$PASTNV),
                         xpre=mean(BASE.DATA$PASTNV),
                         scen=6)

# Scenario 7: Prior Group V Participation
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "Prior Group V Participation", scen=7)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "PASTV",
                         x=mean(BASE.DATA$PASTV) + sd(BASE.DATA$PASTV),
                         xpre=mean(BASE.DATA$PASTV),
                         scen=7)

# Scenario 8: Liberal Democracy Index
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "Electoral Integrity Index", scen=8)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "VDEM_POLYARCHY_LAG",
                         x=mean(BASE.DATA$VDEM_POLYARCHY_LAG) + sd(BASE.DATA$VDEM_POLYARCHY_LAG),
                         xpre=mean(BASE.DATA$VDEM_POLYARCHY_LAG),
                         scen=8)

# Scenario 9: Physical Integrity Index
BASE.FIG1HYP <- cfName(BASE.FIG1HYP, "Physical Integrity Index", scen=9)
BASE.FIG1HYP <- cfChange(BASE.FIG1HYP, "VDEM_PHYSINT_LAG",
                         x=mean(BASE.DATA$VDEM_PHYSINT_LAG) + sd(BASE.DATA$VDEM_PHYSINT_LAG),
                         xpre=mean(BASE.DATA$VDEM_PHYSINT_LAG),
                         scen=9)

# Simulate relative risks with 95\% CI and produce 2 data frame tables, one for Nonviolent outcome the other for Violent outcome
BASE.SIMRR <- mlogitsimrr(BASE.FIG1HYP, BASE.simB, ci=0.95)
RR_NV<-data.frame(rownames(BASE.FIG1HYP$x),BASE.SIMRR$pe[,1], BASE.SIMRR$lower[,1,1], BASE.SIMRR$upper[,1,1])
RR_V<-data.frame(rownames(BASE.FIG1HYP$x),BASE.SIMRR$pe[,2], BASE.SIMRR$lower[,2,1], BASE.SIMRR$upper[,2,1])

##Plot Relative Risks Using Chris Adolph's Tile package
# Set color motif
cols <- brewer.pal(3,"Set1")
#Create tracelines for relative risks for Nonviolent onsets
BASE.TRACERR.NV <- ropeladder(x=BASE.SIMRR$pe[,1],
                                      lower=BASE.SIMRR$lower[,1,1],
                                      upper=BASE.SIMRR$upper[,1,1],
                                      labels=rownames(BASE.FIG1HYP$x),
                                      size=0.65,
                                      col=cols[3],
                                      lex=1.75,
                                      lineend="square",
                                      plot=1
)
#Create tracelines for relative risks for Violent onsets
BASE.TRACERR.V <- ropeladder(x=BASE.SIMRR$pe[,2],
                                   lower=BASE.SIMRR$lower[,2,1],
                                   upper=BASE.SIMRR$upper[,2,1],
                                   labels=rownames(BASE.FIG1HYP$x),
                                   size=0.65,
                                   col=cols[1],
                                   lex=1.75,
                                   lineend="square",
                                   plot=1
)

# Make reference line trace for relative risks (at 1)
vertmarkRR2 <- linesTile(x=c(1,1), y=c(0,1), plot=1)
# Set tick marks for x axis
xat <- c(0.1, 0.2, 0.5, 1, 2, 5)

# Make plot with tile (file format can be adjusted under output/type=)
tile(BASE.TRACERR.NV, BASE.TRACERR.V, vertmarkRR2,
     xaxis=list(log=TRUE, at=xat, labels=paste0(xat,"x")),
     topaxis=list(add=TRUE, log=TRUE, at=xat, labels=paste0(xat,"x")),
     #plottitle=list(labels="Gator 1 compared to (Gator 2)"),
     #xaxistitle=list(labels="relative likelihood of mostly eating \"other\" food"),
     width=list(null=4),
     height=list(xaxistitle=3, plottitle=4),
     gridlines=list(type="xt"),
     output = list(file = "EBCR_Fig1_Ropeladder", family = "Palatino", type= "pdf")
)

###Figure 2: Create Expected Values Simulations for Changes in EPR Status and Groupsize
# Create full factorial set of counterfactuals
groupsize.range <- seq(0,1,by=0.01)          # range of counterfactual groupsizes
SrScen <- JrScen<-ExcScen <- cfMake(BASE.MODEL, data=BASE.DATA, nscen=length(groupsize.range))

#set Groupsize to vary across range under each of 3 different levels of EPR Status (Reference=Senior, Dominant, Monopoly; Junior Partner; Excluded=Self-Excluded, Powerless, Discriminated)
for(i in 1:length(groupsize.range)){
  SrScen <- cfChange(SrScen, "EPR_GROUPSIZE", groupsize.range[i], scen=i)
  SrScen <- cfChange(SrScen, "EPR_STATUS_EXCL", x=0, scen=i)
  SrScen <- cfChange(SrScen, "EPR_STATUS_JP", x=0, scen=i)
  JrScen <- cfChange(JrScen, "EPR_GROUPSIZE", groupsize.range[i], scen=i)
  JrScen <- cfChange(JrScen, "EPR_STATUS_EXCL", x=0, scen=i)
  JrScen <- cfChange(JrScen, "EPR_STATUS_JP", x=1, scen=i)
  ExcScen <- cfChange(ExcScen, "EPR_GROUPSIZE", groupsize.range[i], scen=i)
  ExcScen <- cfChange(ExcScen, "EPR_STATUS_EXCL", x=1, scen=i)
  ExcScen <- cfChange(ExcScen, "EPR_STATUS_JP", x=0, scen=i)
}

#Simulate expected probabilities with 95% CI
SrSimProbs <- mlogitsimev(SrScen, BASE.simB, ci=0.95)
JrSimProbs <- mlogitsimev(JrScen, BASE.simB, ci=0.95)
ExcSimProbs <- mlogitsimev(ExcScen, BASE.simB, ci=0.95)

## Plot expected values for Violent and Nonviolent onsets based on Groupsize and EPR Status across 6 plots

# Create one trace for each predicted category of the response, and each EPR Status
V.Sr.trace <- lineplot(x=groupsize.range,
                       y=SrSimProbs$pe[,2],
                       lower=SrSimProbs$lower[,2,1],
                       upper=SrSimProbs$upper[,2,1],
                       ci=list(mark="shaded"),
                       col=cols[1],
                       plot=4
)

NV.Sr.trace <- lineplot(x=groupsize.range,
                        y=SrSimProbs$pe[,1],
                        lower=SrSimProbs$lower[,1,1],
                        upper=SrSimProbs$upper[,1,1],
                        ci=list(mark="shaded"),
                        col=cols[3],
                        plot=1
)

V.Exc.trace <- lineplot(x=groupsize.range,
                        y=ExcSimProbs$pe[,2],
                        lower=ExcSimProbs$lower[,2,1],
                        upper=ExcSimProbs$upper[,2,1],
                        ci=list(mark="shaded"),
                        col=cols[1],
                        plot=6
)

NV.Exc.trace <- lineplot(x=groupsize.range,
                         y=ExcSimProbs$pe[,1],
                         lower=ExcSimProbs$lower[,1,1],
                         upper=ExcSimProbs$upper[,1,1],
                         ci=list(mark="shaded"),
                         col=cols[3],
                         plot=3
)

V.Jr.trace <- lineplot(x=groupsize.range,
                       y=JrSimProbs$pe[,2],
                       lower=JrSimProbs$lower[,2,1],
                       upper=JrSimProbs$upper[,2,1],
                       ci=list(mark="shaded"),
                       col=cols[1],
                       plot=5
)

NV.Jr.trace <- lineplot(x=groupsize.range,
                        y=JrSimProbs$pe[,1],
                        lower=JrSimProbs$lower[,1,1],
                        upper=JrSimProbs$upper[,1,1],
                        ci=list(mark="shaded"),
                        col=cols[3],
                        plot=2
)

# Plot traces using using Chris Adolph's tile package

tile(V.Sr.trace,
     NV.Sr.trace,
     V.Jr.trace,
     NV.Jr.trace,
     V.Exc.trace,
     NV.Exc.trace,
     #linelabels,
     RxC = c(2,3),
     limits = c(0,1,0,.05),
     #xaxis = list(at=at.x),
     #yaxis = list(at=at.y, major=FALSE),
     xaxistitle = list(labels="Group Size (%)", cex=.8),
     #yaxistitle = list(type="first", labels="Probability of Initiating", x=0.1),
     rowtitle = list(labels=c("Nonviolent", "Violent")),
     columntitle = list(labels=c("Sr. Partner and Above", "Jr. Partner", "Excluded")),
     height=list(columntitle=5),
     width=list(rowtitle=1.5),
     #maintitle = list(labels="Expected Probability of Initiating Campaign", y=1),
     gridlines = list(type="xy"),
     output=list(file="EBCR_Fig2_EVs",width=10, family = "Palatino", type="pdf")
)


###Additional Robustness Checks Included in Appendix#:##

##Producing Table 7: Decade and Regional Controls
#Decade Controls
EBCR$DECADE<-ifelse(EBCR$YEAR<1956, "46-55",
                   ifelse(EBCR$YEAR>=1956 & EBCR$YEAR<1966, "56-65",
                          ifelse(EBCR$YEAR>=1966 & EBCR$YEAR<1976, "66-75", 
                                 ifelse(EBCR$YEAR>=1976 & EBCR$YEAR<1986, "76-85",
                                        ifelse(EBCR$YEAR>=1986 & EBCR$YEAR<1996, "86-95","96-06")))))

DFE.MNM<-multinom(INIT_MN_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + DECADE + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3), data=EBCR, maxit=200)
summary(DFE.MNM)

#Regional Controls
EBCR$REGION<-ifelse(EBCR$GW_COUNTRY_ID<200, "AMERICAS",
                   ifelse(EBCR$GW_COUNTRY_ID>=200 & EBCR$GW_COUNTRY_ID<400, "EUROPE",
                          ifelse(EBCR$GW_COUNTRY_ID>=400 & EBCR$GW_COUNTRY_ID<600, "SSAFRICA", 
                                 ifelse(EBCR$GW_COUNTRY_ID>=600 & EBCR$GW_COUNTRY_ID<700, "MENA",
                                        ifelse(EBCR$GW_COUNTRY_ID>=700, "ASIA","NA")))))

RFE.MNM<-multinom(INIT_MN_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + REGION + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3), data=EBCR, maxit=200)
summary(RFE.MNM)

#Stargazer to generate latex code to produce Table 7
stargazer(DFE.MNM, RFE.MNM,
          title = "Decade and Regional Controls", 
          font.size = "tiny",
          align = TRUE, 
          notes = c("GDP and population lagged 1 year."),
          dep.var.caption = "Odd Models: DV = Nonviolent; Even Models: DV = Violent",
          dep.var.labels.include = TRUE,
          covariate.labels = c(
            "EPR Status: Jr. Partner",
            "EPR Status: Excluded",
            "Group Size",
            "Country Population (log)",
            "Country GDP Per Capita (log)",
            "Prior Group NV Participation",
            "Prior Group V Participation",
            "VDEM Polyarchy Index",
            "VDEM Physical Integrity Index",
            "Decade 56-65",
            "Decade 66-75",
            "Decade 76-85",
            "Decade 86-95",
            "Decade 96-06",
            "Asia",
            "Europe",
            "Middle East and N. Africa",
            "Sub-Saharan Africa",
            "NV YEARS",
            "NV YEARS2",
            "NV YEARS3",
            "V YEARS",
            "V YEARS2",
            "V YEARS3"
          ),
          colnames = TRUE)

##Producing Table 8: Single Equation Models (Nonviolent)

SENV.MODEL<-INIT_NV_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3)

#Logit with and without Clustered SEs
LNV<-glm(SENV.MODEL, family="binomial", data=EBCR)
summary(LNV)
LNV.se<-sqrt(diag(vcov(LNV)))
LNV.gcse <- sqrt(diag(cluster.vcov(LNV, EBCR$EPR_ID)))
coeftest(LNV, cluster.vcov(LNV, EBCR$EPR_ID))
LNV.ccse <- sqrt(diag(cluster.vcov(LNV, EBCR$GW_COUNTRY_ID)))
coeftest(LNV, cluster.vcov(LNV, EBCR$GW_COUNTRY_ID))

#Probit
PNV<-glm(SENV.MODEL, family=binomial(link="probit"), data=EBCR)
summary(PNV)
PNV.se <- sqrt(diag(vcov(PNV)))

#Rare Events Logit using Zelig
SENV.DATA<-extractdata(SENV.MODEL, EBCR, na.rm=T)
REL_ZLG<-zelig(SENV.MODEL, model="relogit", data = SENV.DATA)
summary(from_zelig_model(REL_ZLG))
REL_ZLG.se<-unlist((REL_ZLG$get_se()))

#Stargazer to produce Table 8 in Appendix
stargazer(LNV, LNV, LNV, PNV, from_zelig_model(REL_ZLG),
          se = list(LNV.se, LNV.gcse, LNV.ccse, PNV.se, REL_ZLG.se),
          title = "Single Equations Models: Nonviolent", 
          font.size = "tiny",
          align = TRUE, 
          notes = c("GDP and population lagged 1 year."),
          dep.var.caption = "DV: 1 = Initiate Nonviolent; 0 = No Nonviolent",
          dep.var.labels.include = FALSE,
          #omit = c("NVYEARS", "NVYEARS2", "NVYEARS3", "VYEARS", "VYEARS2", "VYEARS3"),
          covariate.labels = c(
            "EPR Status: Jr. Partner",
            "EPR Status: Excluded",
            "Group Size",
            "Country Population (log)",
            "Country GDP Per Capita (log)",
            "Prior Group NV Participation",
            "Prior Group V Participation",
            "VDEM Polyarchy Index",
            "VDEM Physical Integrity Index",
            "NV YEARS",
            "NV YEARS2",
            "NV YEARS3",
            "V YEARS",
            "V YEARS2",
            "V YEARS3"
          ),
          colnames = TRUE)

##Producing Table 9: Single Equation Models (Violent)

SEV.MODEL<-INIT_V_ONSET ~ EPR_STATUS_JP + EPR_STATUS_EXCL + EPR_GROUPSIZE + POP_LOG_LAG_EXT + GDPPC_LOG_LAG_EXT + PASTNV + PASTV + VDEM_POLYARCHY_LAG + VDEM_PHYSINT_LAG + NVYEARS + I(NVYEARS^2) + I(NVYEARS^3) + VYEARS + I(VYEARS^2) + I(VYEARS^3)

#Logit with and without Clustered SEs
LV<-glm(SEV.MODEL, family="binomial", data=EBCR)
summary(LV)
LV.se<-sqrt(diag(vcov(LV)))
LV.gcse <- sqrt(diag(cluster.vcov(LV, EBCR$EPR_ID)))
coeftest(LV, cluster.vcov(LV, EBCR$EPR_ID))
LV.ccse <- sqrt(diag(cluster.vcov(LV, EBCR$GW_COUNTRY_ID)))
coeftest(LV, cluster.vcov(LV, EBCR$GW_COUNTRY_ID))

#Probit
PV<-glm(SEV.MODEL, family=binomial(link="probit"), data=EBCR)
summary(PV)
PV.se <- sqrt(diag(vcov(PV)))

#Rare Events Logit using Zelig
SEV.DATA<-extractdata(SEV.MODEL, EBCR, na.rm=T)
RELV_ZLG<-zelig(SEV.MODEL, model="relogit", data = SEV.DATA)
summary(from_zelig_model(RELV_ZLG))
RELV_ZLG.se<-unlist((RELV_ZLG$get_se()))

#Stargazer to produce Table 9 in Appendix
stargazer(LV, LV, LV, PV, from_zelig_model(RELV_ZLG),
          se = list(LV.se, LV.gcse, LV.ccse, PV.se, RELV_ZLG.se),
          title = "Single Equation Models: Violent", 
          font.size = "tiny",
          align = TRUE, 
          notes = c("GDP and population lagged 1 year."),
          dep.var.caption = "DV: 1 = Initiate Nonviolent; 0 = No Nonviolent",
          dep.var.labels.include = FALSE,
          covariate.labels = c(
            "EPR Status: Jr. Partner",
            "EPR Status: Excluded",
            "Group Size",
            "Country Population (log)",
            "Country GDP Per Capita (log)",
            "Prior Group NV Participation",
            "Prior Group V Participation",
            "VDEM Polyarchy Index",
            "VDEM Physical Integrity Index",
            "NV YEARS",
            "NV YEARS2",
            "NV YEARS3",
            "V YEARS",
            "V YEARS2",
            "V YEARS3"
          ),
          colnames = TRUE)
